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Maxwell equations are solved in a layer comprising a finite number of homogeneous isotropic 
dielectric regions ended by anisotropic perfectly matched layers (PMLs). The boundary- value 

>v 

43 ' problem is solved and the dispersion relation inside the PML is derived. The general expression 

of the eigenvalues equation for an arbitrary number of regions in each layer is obtained, and 

(N . 

both polarization modes are considered. The modal functions of a single layer ended by PMLs 

IT) 

£\J ' are found, and their orthogonality relation is derived. The present method is useful to simulate 

o: 

scattering problems from dielectric objects as well as propagation in planar slab waveguides. 
Its potential to deal with more complex problems such as the scattering from an object with 

O ' 

""^5 ' arbitrary cross section in open space using the multilayer modal method is briefly discussed. 
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1 Introduction 



The perfectly matched layer (PML) is a fictitious material which does not reflect incident prop- 
agating waves regardless of the incident angle, frequency and polarization. It was firstly intro- 
duced by Berenger ^ |2] as a useful absorbing boundary condition to truncate the computational 
domain in finite-difference time-domain applications. A different formulation of the PML, given 
by Sacks et al. [3], is based on exploiting constitutive characteristics of anisotropic materials to 
provide a reflectionless interface. This formulation offers the special advantage that it does not 
require modification of Maxwell equations [3] . Both formulations of the PML are very popular 
among the electrical engineering community and their use in the optics community has been 
growing in the last few years H|-|H|. PMLs look attractive, mainly because of their potential 
capacity to simulate open-structure problems by bounded computational domains, which conse- 
quently reduce significantly the computational costs involved in the calculations Even 
though this kind of medium is usually used in the framework of finite-element methods |12|-|14j. 
there are also studies which incorporate the anisotropic absorber in different approaches [7j, |15j - 
|19| . In particular, the modal approach appears as an interesting alternative for the description 
of the fields in index-modulated structures since it highlights the physics of the problem. 
The modal approach has been applied by many authors to dielectric lamellar gratings in classical 
[^nj-jH] as well as in conical mounting |231 I24j . In his work, Li [21] derived rigorously the 
eigenfunctions and also their completeness and ortogonality relations from the boundary-value 
problem. Later on, this formalism was extended to deal with multiply grooved lamellar gratings 
|25| . where the eigenvalues equation for an arbitrary number of different dielectric zones was 
obtained and carefully analyzed for real and complex refraction indices. The modal formulation 
presented in |25| was then extended to deal with finite structures such as index-modulated 
apertures 

A further step on the development of the modal method was the application of the multilayer 
approximation [2^ not only to infinite gratings j^Tj but also to finite structures with arbitrary 
shapes of the corrugations [21]- [30]. However, all these works were restricted to structures 
laterally closed by perfect conductors. If illuminated by finite beams, these structures can 
approximately simulate the scattering from purely dielectric finite structures |28|-|3U|. 
It is well known that either the pseudeperiodic condition in infinite gratings or the assumption 
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Figure 1: (a) Scheme of the scattering problem by an index-modulated aperture; (b) Scheme of 
the propagation problem in a planar slab waveguide. 

of perfect conductivity in finite structures forces the eigenvalues set to be discrete. On the 
other hand, open problems like propagation in dielectric waveguides or scattering from dielectric 
objects in open space, have a continuum set of eigenvalues. Therefore, it is airways convenient to 
limit the problem domain, in order to avoid large volumes of calculus. This can be done by means 
of perfectly matched layers, which ensure absorption and attenuation of the incident radiation 
and consequently simulate the open space better. A detailed formulation of the discretization 
of the continuous spectrum by means of a Dirichlet boundary condition for the complex plane 
has been recently presented where the solution of the Sturm-Liouville problem for a layered 
structure is based on the PML boundary condition. The completeness of the eigenmodes of 
parallel plate waveguides with PML terminations has also been proven |32| . 
The purpose of this paper is to find the analytical expressions of the eigenmodes of a single index- 
modulated layer ended by PMLs at both sides, and also to derive the orthogonality relations 
that satisfy the set of eigenfunctions. A modal formulation is presented to apply the eigenmodes 
expressions to the solution of the scattering problem by a two-dimensional object of rectangular 
cross section in open space. We cover simultaneously the eigenmodes of an index-modulated 
aperture (schematized in Fig. la) as well as those of a planar slab waveguide (Fig. lb), since 
both problems have an equivalent formulation. Numerical examples comparing the resulting 
eigenvalues with those found by a different approach are shown. The present paper is more 
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general than Ref. ^7] since we find the eigenvalues and the eigenmodes of a layer comprising 
an arbitrary number of different homogeneous regions, with the aim of extending the modal 
method to more general structures ended by PMLs. 

Finally, a brief discussion on the potential applicability of the modal method to the scattering 
problem by a 2-D object of arbitrary cross section in open space is given. This extension can 
be made by means of the multilayer approximation |2J , and the application of any propagation 
algorithm such as the R- matrix pi 171 127] [35] -[55]. 

2 Modal formulation 

In this section we develop a rigorous method to find the exact modes of a structure as that 
of Fig. 2. The modes of such a structure can be then used to find the solution of waveguide 
problems and scattering problems. 
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Figure 2: Configuration of the problem. 

Consider an index-modulated layer invariant in the z direction, as shown in Fig. 2. Each 
homogeneous zone Xj < x < xj + Aj of width A,- has a permittivity €j and a permeability fij, 
where e,- and fj,j can be complex numbers. The layer is ended at both sides by perfectly matched 



4 



layers backed by perfect conductors. The PML is a particular anisotropic medium especially 
designed to absorb the received radiation [5J , and its constitutive parameters are epML = eA and 
Apml = where e and [x are the permittivity and permeability of the adjacent isotropic zone, 
respectively. For an interface parallel to the (y, z)-plane, the tensor A must be defined as: 



A 



1/6 
6 
6 



(1) 



The parameter 6 is a complex number: its real part is related to the wavelength inside the 
PML, and the imaginary part accounts for the losses in the material, i.e., the attenuation of the 
propagating waves. 

As stated above, to find the eigenmodes of the structure depicted in Fig. 2 can be useful for 
waveguide and scattering problems. We are interested in studying the propagation of plane 
waves along the structure, and then a time dependence e~ %wt (uj being the frequency of the 
plane wave) is implied and suppressed in all the paper. Since the structure is invariant in the 
z-direction, the problem can be separated into two scalar problems for E z (TM polarization) 
and H z (TE polarization). Then, the x- and y- components of the fields can be written in terms 
of the z-components [55] . 
Within the PML, Maxwell curl equations are 



V x E = i-uAH , 

c 



(2) 



V x H 



-i-eAE , 

c 



(3) 



where c is the speed of light in vacuum and E and H are the electric and magnetic fields, 
respectively. In the case of a slab waveguide (Fig. lb), we are interested in finding the solu- 
tions of Maxwell equations that have a z-dependence of the form e tkzZ . Then, combining eqs. 
© and (J5J), and taking into account the invariance of the problem along the z direction, we 
get the propagation equations for a layer with x-dependent constitutive parameters, for both 
polarization modes: 



fi(x) d 
b(x) dx 



1 



dE z (x,y) 



b(x)fj,(x) 



dx 



d 2 E z (x,y) 9 . . 

|A_^Z +7 2 (x)^(x,y) = , 



(4) 



e(x) d 
b(x) dx 



1 

b(x)e(x) 



dH z (x,y) 
dx 



+ 



d 2 H z (x,y) 
dy 2 



+ -f 2 (x)H z (x,y)=0 



(5) 
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where 



e{x)ji(x) — k 2 , 



(6) 



and b(x) = 1 in the isotropic regions. For the scattering problem of a plane wave with wave 
vector k in the (x,y)-plane impinging on the structure (Fig. la), k z = and then 7 2 (x) = 
k 2 (x) = (u> 2 /c 2 )e(x)fi(x). On the other hand, for a slab waveguide the problem is invariant in 
the y direction, and then d y = 0. Since the form of eqs. @ and © is identical, we unify the 
treatment of both polarizations in a single differential equation 



<j{x) d 
b(x) dx 



1 



dip(x,y) 



b(x)a(x) dx 



o 



where 



fi(x) for TM polarization 



(Tlx) 



(7) 



(8) 



e(x) for TE polarization 
and ip represents either E z (TM case) or H z (TE case). 

The eigenmodes of the structure are the set of linearly independent functions that satisfy by 
themselves the boundary conditions at all the interfaces x = Xj, j = 0, ...,N + 1, and form a 
complete basis. In particular, since the layer is ended by a perfect conductor, we require that 
the tangential component of the electric field vanishes at x = xq and at x = L. To do so, we first 
solve eq. ® in each region (isotropic and anisotropic) and then match these partial solutions 
correspondingly. For the most general case in which d y ^ a separated solution is proposed: 



^(x,y) = X(x).Y{y) , (9) 

and substituting eq. @ in eq. (JJJ) we get two ordinary differential equations for the functions 
X(x) and Y(y) 

d 2 Y 



dy* 



+ v 2 Y 



, 



(10) 



a(x) d 
b(x) dx 



1 



dX 



b(x)a(x) dx 



+ ( 7 2 (x) -v 2 )X 



, 



(11) 



where v is a constant. The formal solution of eq. (|10|) is straightforward 



Y(y) = acos(vy) + bsin(vy) 



(12) 
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and we will now focus on the differential equation for X{x). Notice that for the scattering 
problem j 2 (x) = k 2 (x) — v 2 and for the slab waveguide problem 7 2 (x) = k 2 (x) — k 2 . Then, we 
can unify both cases by setting p = v 2 in the first case and p = k 2 in the second case. 
Eq. ifTTj) together with the boundary conditions at the ends of the layer (x = x$ and x = L) 
pose a boundary-value problem 

CX = pX , (13) 



where C is the differential operator 

C 



a(x) d \ 1 d 1 lo 2 . . , . 

+ —p{x)e(x) . 14 
cr 



b(x) dx [b(x)a(x) dx 

Since b is a complex number, the operator £ is a non-self-adjoint operator [SZj. Prom the theory 
of non self-adjoint problems we know that the eigenvalues are infinite and complex in general, 
and the eigenfunctions do not necessarily form a complete and orthogonal set [H7j. To find a 
useful solution to our problem, is then necessary to consider the adjoint problem of (|T3|) 

C*X^ = ptxt , (15) 

where the asterisk denotes complex conjugate and the ' denotes adjoint. The eigenvalues of the 
adjoint problem are the complex conjugates of those of the direct problem, i.e., 

P ] = P* , (16) 

and the sets of eigenfunctions {X m (x)} and {xt(x)} form a bi-orthonormal set such that 
(X m ,Xn) = S mn , where the internal product (X m ,Xn) is defined as 

(X m ,Xl) = [ ^l Xm (x)(Xi(x)Tdx . (17) 
J <?{x) 

Then, any continuous and piecewise differentiable function that satisfies the boundary conditions 
at the ends of the layer can be expanded as 

f(x)=J2 A mX m {x) . (18) 
m 

Taking into account these facts, we are going to find the solutions of eq. Qll|h In each homoge- 
neous region (j), we propose a solution of the form 

Xj(x) = Aj cos[uj(x — Xj)] + Bj s'm[uj(x — Xj)] for Xj < x < Xj + Aj , (19) 
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and substitute it into eq. to get the dispersion relation 



u 2 

t| + P = 7|, (20) 
j 



for the two kinds of zones we have in the structure: isotropic (bj = 1) and anisotropic (bj £ C, 
bj ^ 1). Notice that p is a constant for each eigenmode, and the constitutive parameters of 
the PML satisfy eo = eiA, po = Pi^-, and then = k\. The same occurs at the other end of 
the layer, and therefore, uq = bu\ and um = bujy-i, i.e., the eigenvalue in the PML region is 
b times the eigenvalue in the adjacent isotropic region. Then, the solution of eq. in the 

whole layer can be expressed in terms of harmonic functions, where the amplitudes are such that 
satisfy the boundary conditions at the vertical interfaces, and the eigenvalues u and p satisfy 
the corresponding dispersion relation, depending on the zone (eq. (f2T))0 . Defining X(x) as: 

~ IX ' 

X(x) = — [—Bj cos [uj(x — Xj)] + Aj sin [uj(x — Xj)]] for Xj < x < Xj + Aj , (21) 
it is easy to verify that |26j : 



a 



X(x) = X(xj) cos(uj(x — Xj)) -X(xj) sin(uj(x — Xj)) , for Xj < x < Xj + Aj 

Uj 

~ — IX ' 

X(x) = X(xj) cos(uj(x — Xj)) H — -X(xj) sm(uj(x — xj)) , for Xj < x < Xj + Aj (22) 

CTj 

Notice that X(x) is proportional to the normal derivative of X(x) at the vertical interfaces, 
which implies that it represents the tangential component of the magnetic (electric) field in the 
case of TM (TE) polarization. 

Equations (|22jl provide us with a relation between two field quantities at both sides of each 
homogeneous zone bounded by the interfaces at x = Xj and x = Xj + Aj. This relation can be 
expressed in matrix form as: 

X(x j+1 ) \ / X(xj) \ 

= Mj(p) , (23) 

X(x j+1 ) J \ X(xj) J 

where M,-(p) is a matrix given by: 

cos (uj Axj) — sin(ujAxj ) aj / Uj 
sin(ujAxj) Ujjoj cos(tijAxj) 

X(x) and X(x) represent vectors containing the modal functions X m {x) and X m (x), respectively, 
and u and p are related by the dispersion relation. The subscript m denoting the mode has been 



Mj(p) 



(24) 



suppressed as it is understood that relations Q22 |) -(|24 |) hold for each one of the modal terms. 
This procedure can be applied N times to get a relation between the fields at the perfectly 
conducting walls at x = xo and x = xl- In such a case we have a relation of the form: 

/ X(x L ) \ / X(x ) \ 

= M(P) , (25) 

V X(x L ) J \ X(x ) / 

where M(p) is a product matrix: 

M(p) = M L (p) M L _i(p) . . . M 2 (p) Mi(p) . (26) 
M(p) is a 2 x 2 block matrix. This is a well known result from the theory of stratified media 

Imposition of the boundary condition at xq and at xl, i.e., the tangential component of the 
electric field must vanish on the surface, yields a condition on the non-diagonal block elements: 

M 12 (p) = for TM polarization (27) 
M 21 (p) = for TE polarization . (28) 

These conditions determine eigenvalues equations for p, that must be solved by means of numer- 
ical techniques. These equations have already been studied for structures formed by isotropic 
regions in the case of an infinite periodic grating and of an index-modulated aperture |261l29j . 
The expressions of the equations found in this work reduce to those in Refs. |261 129j for 6 = 1. 
Also, a recursive formula to find the eigenvalues equation of a layer with N different regions could 
be derived in the present case [25]. Finally, the eigenmodes of the layer bounded by perfectly 
matched layers backed by perfectly conducting walls are given by 

X m (x) = A m jCos[u m j(x - Xj)] + B m jsm[u m j(x - xj)] for Xj < x < xj + Aj , (29) 

where u m j satisfy the corresponding eigenvalues equation (|2*7j) or (|28|) . The modal coefficients 
A m j and B m j are easily obtained by successive application of eqs. (|22j) 

A m j = A m>j -i cos[u m j-i(xj - Xj-i)) + -B m ,j-i sm[u m j-i(xj - Zj-i)] 

Bmj = _!Mzi_^L_ {B m>J -_i cos[n mj „i(xj - Xj-i)] - An,j-i sm[u m ,j-i(xj - Xj-i)]} (30) 

Umj &j — 1 
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3 Examples 

The eigenvalues equation is the critical part of any modal formulation, particularly for layers 
comprising different kinds of materials. Equations of this kind have been solved in Refs. |25| . 
|26| and [2S| f° r rea l refractive indices (non lossy dielectrics) and in |30| for lossy materials. 
A detailed study of the respective equations has been performed in the above works to avoid 
numerical problems with the rootfinding algorithms. In |26|.|28| and |29| the authors use the 
separation constant along the y direction (v 2 ) as the variable of the equation, whereas for highly 
conducting metals, i.e. regions of refraction index with a large imaginary part, it is convenient 
to use the separation constant in the x direction (u) as the variable .30]. This behaviour can 
be understood taking into account the location of the eigenvalues in two limit cases: (i) ideal 
dielectric isotropic material (real refraction index), and (ii) perfect conductor (infinite refraction 
index). In the first case, the eigenvalues (v 2 ) are real, and then, for small perturbations of 
the refraction index (lossy dielectrics) a small departure from these values is expected. On the 
contrary, for the perfect conductor the eigenvalues (u) are real and can be calculated analytically 
|39| . Therefore, we expect the eigenvalues corresponding to a highly conducting metal not to 
depart too much from those real values. 

As a first example, we consider the case of a symmetric structure formed by three zones (PML 
- dielectric - PML) of widths Ap and Ai, and then u$ = U2 = up. This example actually 
simulates the open space, and can be useful to study the coupling between a slab waveguide and 
open space |17j . For this case, the direct application of the procedure described in the previous 
section yields an eigenvalues equation for both polarization modes, that after some manipulation 
can be reduced to: 

sin(2« P Ap + u l Ai) = . (31) 
Taking into account that up = u\b, the solutions of eq. 1)31(1 are 

«im = 77T^ 1 — a~~ ' for m mt eger. (32) 

Z LAp + L\\ 

The eigenvalues obtained in (|32|) are the same as those found by Derudder et al. by means of a 
coordinates transformation in the propagation equation (eqs. (16) and (17) in Ref. |17j). 
It is important to remark that an increase in the number of homogeneous regions in the layer 
produces a more complex eigenvalues equation, since it is obtained from the determinant of a 
matrix product (eq. (|2fi)l). The number of terms and the number of trigonometric factors in 
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Figure 3: Comparison of the eigenvalues obtained by the present method and that of Derudder 
et. al (eqs. (24) and (25) in Ref. [T7]) for a structure with 5 zones: A P /X = 0.04, Ai/A = 
A 3 /A = 0.5, A 2 /A = 2, ni = n 3 = 1, n 2 = 1.3, 6 = 2 — 2*, TM polarization. 

each term increases, causing abrupt oscillations in the equation. However, since we are focusing 
on the solution of scattering problems from homogeneous objects in open space, the number 
of different regions to be considered in this case, including the PML, is restricted to five. The 
explicit expressions of the eigenvalues equations for five zones and for Ai = A3 and jij = /ioY? 
are given in the Appendix for TE and TM polarizations. As it can be observed in Fig. 3 for 
TM polarization, the solutions of the equation are coincident with those obtained from eqs. (24) 
and (25) in Ref. ^7]. These equations correspond to the eigenvalues equations for odd and even 
modes, respectively. 

In Fig. 4 we show the evolution of the eigenvalues as the imaginary part of the PML parameter 
b is increased, for a structure comprising five zones and TE polarization. For Im(b) = 0, the 
eigenvalues u\ are real, and their values depart from the real axis as Im(b) is increased. 
As it has been stated above, any system of planar regions can be studied by the modal method. 
If we take into account that any profile can be approximated by a stack of rectangular lay- 
ers (multilayer approximation) the scope of applicability of the modal method broadens 

11 



1 - 



- 



-1 - 



-2 - 



D □ 








• 

o • 


□ 


□ 


□ 


o 


• 






o 


• 




• 



E "3- 
-4 - 
-5 - 
-6 - 



■ 


lm(b)=0 


□ 


lm(b)=0.05 


• 


lm(b)=0.1 


o 


lm(b)=0.15 



— r 

20 







40 



60 



80 



100 



Re( Ul ) 



Figure 4: Evolution of the eigenvalues with the imaginary part of the PML parameter b, for a 
structure with 5 zones: Ap/A = Ai/A = A2/A = A3/A = 0.2, ri\ = 723 = 1, n<i = 1.5, Re(b) = 1 
and TE polarization. 

significantly. Since we have already found the exact expression of the field in terms of the eigen- 
modes of a rectangular layer, the method can be applied to a multilayer system by means of any 
propagation algorithm such as the T-matrix, R- matrix, and the S-matrix |27| . The procedure 
is essentially that formulated in Refs. |28|.|29| for dielectric apertures. The main difference is 
that in the present case, i.e., when the layers are ended by PMLs, the definition of the modal 
functions U^^x) (eq. (9) in Ref. j22]) should be substituted by the expression of X m (x) given 
in eq. (|29j) for each layer. 

4 Summary and conclusions 

Analytical expressions for the eigenmodes of an inhomogeneous layer ended by particular aniso- 
tropic regions (PMLs) have been found. The general method was described for an arbitrary 
number of regions of different widths and dielectric materials in the structure. The eigenvalues 
equation for structures comprising lateral anisotropic lossy regions was found, what opens up the 
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possibility of application of the modal method to simulate scattering problems from 2-D objects 
in open space. A complete set of modal functions was obtained, and also the orthonormality 
relation between them was derived. Simple examples of eigenvalues equations for planar slab 
waveguides (three and five zones in the layer) have been presented, and the results have been 
compared with those found in the literature. The present formulation constitutes an extension 
of the modal method and it is a basic tool that can be applied to more general problems such 
as propagation in waveguides of arbitrary cross section and scattering problems from arbitrarily 
shaped objects. 
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Appendix 

Explicit expressions for the eigenvalues equations obtained for a symmetric structure of five 
zones (PML - 1 - 2 - 1 - PML), with fij = (i \/j. 
TM polarization 



2cos(itpAp) cos(wiAi) sin(/upAp)[2sin(«2A2) 2 — 1] + 
2sin(n2A2) cos(ti2A2) cos(itiAi)[2sin(«pAp) 2 — 1] + 
2 cos(upAp) sin(u2A2) cos(u2A2) sin(uiAi) sin(iipAp) ( — H — - ) + 




sin(«pAp) 2 sin(uiAi) 




cos(upAp) 2 sin(uiAi) — cos(u2A2) 2 — + sin(ti2A2) 2 — =0 



(33) 
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TE polarization 



2sin(iipAp) cos(upAp) cos(uiAi)[l — 2sin(u2A2) 2 ] + 
2sin(it2A2) cos(«2A2) cos(uiAi)[l — 2 sin(itpAp) 2 ] — 

/ U\ 62 V>2 ^1 \ 

2cos(upAp) sin(n2A2) cos(n2A2) sin(niAi) sin(npAp) 1 

V. ''2 ' : '': (-2 J 

U2 £l sin(niAi) [sin(npAp) 2 cos(u2A2) 2 + sin(n2A2) 2 cos(npAp) 2 ] + 
1 62 sin(uiAi) [cos(n2A2) 2 cos(npAp) 2 + sin(u2A2) 2 sin(upAp) 2 ] = 
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